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We report the results of three-dimensional molecular dynamics simulations of sheared granular 
system in a Couette geometry^ The simulations use realistic boundary conditions that may be 
expected in physical experiments. For a range of boundary properties we report velocity and density 
profiles, as well as forces on the boundaries. In particular, we find that the results for the velocity 
profiles throughout the shearing cell depend strongly on the interaction of the system particles 
with the physical boundaries. Even frictional boundaries can allow for significant slippage of the 
particles, therefore, reducing the shear in the system. Next, we present shear rate dependence of 
stress, including mean force and force fluctuations, both for controlled volume, and for controlled 
stress configurations. We discuss the dependence of solid volume fraction on shear rate under the 
constant pressure condition, and Bagnold scaling in volume controlled experiments. In addition, we 
study the influence of oscillatory driving on the system properties. 



I. INTRODUCTION 

A significant understanding of granular flow results 
from experimental measurements and numerical simula- 
tions. Simulations have proved to effectively complement 
the experimental measurements. Among many known 
advantages of computer simulations in any field of physics 
we can distinguish the following three types: 1) The abil- 
ity to study the regions of parameter phase space that are 
difficult to access in experiment. For example, regarding 
granular experiments, high-frequency and/or high am- 
plitude driving can be easily simulated but is difficult to 
achieve in experiments because of considerable power re- 
quired; 2) The ability to study the system under (often 
simplified) conditions that are not reachable or are very 
difficult to have in experiments. For example, 2D simu- 
lations as opposed to 3D experiments, periodic boundary 
conditions as opposed to solid walls boundaries, etc.; 3) 
The ability to obtain more complete and detailed infor- 
mation under the same or close to the same conditions 
as in an experimental study. 

However, in the subfield of sheared granular flow, most 
of the known numerical studies exploit either the first 
or the second advantage of simulations from the above 
ij g ^ | 2 J 3 J 4 J 5 J 6 J 7 J 8 J 9 j n p ar ticular, there are very few efforts 
to model realistically the physical boundaries, and this 
is the most common discrepancy between the set-up of 
laboratory experiments and the set-up of numerical sim- 
ulations. This effect of physical boundaries, such as the 
walls of the container, is often and rightfully considered 
responsible for the discrepancy between the results of the 
experiments and of the simulations, or between the pre- 
dictions of the existing theories and the results of ex- 
periments or simulationsii2iii Jenkins and Richman», 
for example, calculated boundary conditions in a spe- 
cific limit for plane flow of identical, smooth, inelastic 



disks interacting with a bumpy wall. Louge, Jenkins and 
Hopkins^ and later Lougei^ tested these theoretical pre- 
dictions for rapid sheared granular flows using computer 
simulations. The work of Nott, Alam, Agrawal, Jackson 
and Sundaresanei ali^ presents the theoretical study of 
the effect of boundaries on the plane Couette flow, in- 
dicating the possibility of many different stable and un- 
stable states of the flow, completely determined by the 
properties of the boundaries. However, to the best of our 
knowledge, none of the simulations reported in the liter- 
ature for sheared granular flow are using the boundary 
models that closely reflect the geometry and the con- 
ditions of a Couette cell. Instead, simple models and 
boundary conditions are conventionally used to set and 
study sheared granular flow. Campbell and Brennen-, 
for example, recognize the effect of shearing wall proper- 
ties on velocity and granular temperature profile in 2D 
shearing cell by presenting in parallel the results for two 
models of particle- wall collisions: Type A model assumes 
the interaction between particles and the wall is similar 
to the interaction between particles and particles, and; 
Type B model attempts to approximate a no-slip condi- 
tion by setting the after-rebounce velocity of a colliding 
particle equal to the velocity of the wall. Type B model 
has then been studied more extensively by many authors 
because it provides a constant shear rate and constant 
volume fraction throughout the sample and thus makes 
a convenient test case when comparison with theoreti- 
cal predictions is considered. However, the experimental 
data indicate strongly the presence of boundary effects 
on the shear flow. The examples are the formation of a 
shear band, or other non-linear shear rate and volume 
fraction distributions in the sample (see e.g. Ref . ITEf) . or 
the drastic difference between the shearing using a wall 
with glued particles on it as opposed to shearing the sys- 
tem by a plane frictional wall. Similarly, the periodic 
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boundary conditions, that are often used in simulations, 
neglect the effects of stationary walls on the shear flow. 

In this paper we report the systematic study of the 
effect of boundary conditions on sheared granular flow 
in a Couette cell in zero gravity using event-driven sim- 
ulations. The paper is organized as follows. Section ITTI 
describes the numerical approach and relevant parame- 
ters of simulations. In Sec. II I II we present the results 
for bulk velocity and volume fraction profiles for the va- 
riety of boundary conditions, such as the system with 
rough shearing wall with and without glued particles on 
it. Also we consider simulations with oscillating bottom 
wall. We show the results for the systems of varied sizes, 
varied properties of the walls of a Couette cell, varied in- 
tensity of external driving, and discuss the effect of each 
varied parameter on the profiles and equilibrating times. 
In Sec.|W]we discuss stresses on the physical boundaries 
of a Couette cell and their distributions. We compare 
the normal and shear component of the stress distribu- 
tions and discuss the factors determining the widths of 
the distributions and their average values. 

All the results in Sees. IIII AlIVl are obtained for the 
systems of fixed or directly controlled total volume. In 
Sec. we present the results for the stresses, velocities 
and volume fraction profiles for the systems with stress- 
controlled boundaries. We then discuss the differences 
between volume controlled and stress-controlled simula- 
tions. In particular, we consider the different response to 
the increase of the shearing velocity in these two config- 
urations. 



II. 



SIMULATION DETAILS 



Numerical algorithm of our choice, event-driven algo- 
rithm for hard spheres, is well described elsewhere, for 
example in Ref. Il6l In zero gravity all particles follow 
linear trajectories between the collisions. All collisions 
between particles are instantaneous and binary. Consider 
two colliding particles with diameter <j\ and 02, masses 
mi and 1712, positions ri and T2, linear and angular ve- 
locities vi, oj\ and V2, u>2- Then velocities after collision 
are: 
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Here n = (Y2 — i"i)/|r2 — ri| is the normal unit vector, v„ 
and v t are the normal and tangential components of the 
relative velocity v c = vi — v 2 — [(^cti uj\ + \a2UJ2) x n] 
of particles at the contact point. Total linear and angu- 
lar momentum are conserved during a collision, however, 
total translational and rotational energies are lost. 



Energy dissipation is controlled by three parameters^ 
the coefficient of restitution e, coefficient of friction /j, 
and coefficient of tangential restitution (3. The algorithm 
is adapted to include the particle's interaction with phys- 
ical boundaries in the following way. The collisions be- 
tween particles and the top or bottom walls (also called 
horizontal walls) assume that the wall is a particle of in- 
finite radius, infinite mass, and of linear velocity equal 
to the wall's velocity at the point of contact. Three dis- 
sipation parameters, e w fi w and (3 W , characterize these 
collisions. In the collisions between free and glued (to be 
explained below) particles we use the same dissipation 
parameters as in the collisions between free particles. In 
the calculation of the re-bounce velocity we assume that 
a participating glued particle has an infinite mass. Fi- 
nally, the side walls are the particles of infinite mass and 
of infinite radius and with the surface normal at the point 
of contact pointing horizontally toward the center of the 
cell. The properties of the side- walls are characterized by 
the separate set of dissipation parameters, e s , /j, s and j3 s . 

Experiments and theoretical studies show that the co- 
efficient of restitution e noticeably depends on the impact 
velocityi 18 i 19 To account for this effect, we set this coef- 
ficient for all types of collisions to be velocity dependent 
in a way suggested in Ref. I20I 
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Here v n is the component of relative velocity along the 
line joining particle centers, B = (1 — e)vy 3 ^ 4 \ vq ~ 
100 (d) I sec e is a restitution parameter that characterizes 
the material properties of the particles or walls: Decreas- 
ing e is equivalent to increasing the energy dissipated 
during each collision. 

Coefficient of tangential restitution, (3, defined as the 
ratio of the tangential components of the relative veloc- 
ity after and before collision, is assumed to be the one 
suggested in Ref. \2 ll and used in Ref. l20l i.e. 



— 1 + + e)(7/2)v n /v t , for sliding contacts 
0o, for rolling contacts 
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Here /j, is the coefficient of friction, e is given by and 
/3o is a parameter which also has a meaning of restitution 
of velocity in the tangential direction for sliding contacts. 
As in Ref. 2jJ, the choice between rolling and sliding so- 
lution depends on the ratio of normal to tangential com- 
ponent of relative velocity, and on the parameters flo, e, 
(i. The same functional dependence is assumed for (3 W , 
e w , and for /3„, e s . 

Table [I] summarizes our choice of dissipation parame- 
ters for typical configurations. When different dissipa- 
tion parameters are used in the simulations, the results 
are always compared to the typical configuration. 

We have carried out the simulations using a 3D gran- 
ular system of polydisperse inelastic frictional spheres in 
the Couette cell, Fig. ^ with stationary side walls, ro- 
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TABLE I: Typical choice of dissipation parameters. 



particles (including glued) 


e = 0.6 


H = 0.5 
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frictional inelastic particles 
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very frictional and inelastic 
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Pos = 0.35 


slightly frictional and inelastic 



y axes 




Top shearing wall 



Bottom wall 




FIG. 1: Granular matter between two concentric stationary cylinders - side walls (not shown in the figure). Top wall is rotating 
around the vertical y axes. Bottom wall can be vibrated in the vertical direction in volume controlled simulations, or it can 
move in the vertical direction in stress controlled simulations. 



tated top wall, and the bottom wall that can move in the 
vertical direction. 

The last feature allows us to oscillate the bottom wall. 
The motivation for these oscillations is that they are often 
used in laboratory experiments performed with granular 
systems of high volume fraction in order to fluidize other- 
wise jammed stated. We use the oscillating bottom wall 
to study the effect of oscillations on the granular flow and 
to make the connection with the laboratory experiments, 
where possible. 

In addition, the moving bottom wall feature is ac- 
tively used in stress controlled simulations, where we al- 
low granular particles to determine the position and the 
velocity of the bottom wall on their own. Here the stress 
constraint determines the average volume fraction. We 
refer to Section for more detailed discussion of this 
configuration. 

The particles are polydisperse with diameter values 
randomly distributed in the range [0.9,1.1] (d), where 



(d) is the average diameter of particles, that is used as 
a natural length scale. Our typical configuration consist 
of N — 2000 particles, and fixed radii of inner cylinder 
and outer cylinder (Ri = 8(d), R = 12(d)). The height 
H of the Couette cell is variable and in general is a func- 
tion of time and preset average total volume fraction v 
for the volume controlled simulations, and it is the func- 
tion of time and preset average stress on the bottom wall 
in stress controlled simulations. In the case of volume 
controlled simulations we have 

H{v,t) = ^4.22- +Asin(ujt)j (d). (6) 

Here A, u> are the amplitude scaled by (d), and the 
frequency in rad/sec of the bottom wall vibrations, if 
any. Thus, in absence of vibrations, H = 10.54(d) for 
v = 40%; H = 21.10(d) for v = 20% and H = 8.44(d) for 
v = 50%. Our typical choices for the amplitude and fre- 
quency of the vibrations of the bottom wall are A = 1(d) 
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and uit — 230 rad/ sec or / = w t /27r = 36.6 Hz. These 
amplitude and frequency set the dimensionlcss parame- 
ter r = Aui 2 / g to the value which is above critical for the 
onset of fluidization, r = 5 > T c = 1. We use the value of 
gravitational constant g = 9810(d) /sec 2 (= 981 cm/ sec 2 
when (d) =0.1 cm). 

Another system size studied in this work is a wide con- 
figuration that has N = 5000 particles in a Couette cell 
with radii of inner and outer cylinders Ri = 5(d) and 
R Q = 15(d), respectively, and H given by ©■ 

We assume the surface properties of side walls to be 
the same for inner and outer cylinder and different from 
the properties of top and bottom walls. The top wall 
can be covered with particles glued to its surface. These 
glued particles have the same polydispersity, mean diam- 
eter and material properties as the free particles. They 
are glued into the dimples on the surface so that each 
one protrudes inside the cell by the distance of half it's 
diameter. Figure [3 shows a typical distribution of glued 
particles on the inside surface of the top wall. For both 
typical and wide configurations the number of glued par- 
ticles corresponds to « 30% surface fraction. 

III. VELOCITY AND VOLUME FRACTION 
PROFILES 

To measure velocity and volume fraction profiles, the 
volume of Couette cell is divided in N s ring slices, Fig. [3| 
(a). Each slice is assigned a y bin number i y = 1, N s . 
In most cases, we set N s = 10 and the height of one 
ring slice is h y — 1.05(d). When the bottom wall is os- 
cillated, we account for variable volume of the bottom 
slice or slices. In our measurements we average a studied 
quantity both over the volume of a y bin and over the 
time interval t p , typically t p = 0.5 sec. Within this time 
interval we measure the particle's velocities and volume 
fractions periodically every St p leading to N p — t p /St p 
different distributions. It is convenient to express times 
in terms of the following time scale 



Time T w correspond to the period of oscillation with fre- 
quency uj t . Thus t p — 18.3 T^,, and we use St p = T w /10. 
T w is used as a scale independently of whether the sys- 
tem is vibrated or not. All particles are binned in the 
corresponding ring slices, and the average solid volume 
fraction and the average tangential velocity in each bin 
are calculated as a function of y bin number. 

For wide configurations, in addition we measure the 
volume fraction as a function of the distance to the inner 
cylinder. In this case the volume of Couette cell is divided 
into a set of vertical shells as in Fig. |3Jb). Each shell 
is assigned r bin number. The averaging procedure is 
analogous to the one described above for y profiles. 

In the figures of this paper, the steady state veloc- 
ity profiles are shown as lines with filled triangles. The 



"shadow of broken lines" below the steady state veloc- 
ity profile are transient state velocity profiles measured 
periodically - typically every 1/12 of the total simula- 
tion time. Filled circles mark the steady state volume 
fraction profiles. Transient states for volume fraction are 
also shown as the collection of broken lines. For each 
case we show the scaled energy plot as a function of time. 
The scaled energy is measured kinetic energy per particle 
scaled by the energy of a particle of average diameter and 
of average linear top wall velocity. These energy plots al- 
low to determine whether and when the system reaches 
a steady state. 



A. Typical configurations with inelastic and 
frictional sidewalls 

Figure 0] shows the velocity profiles, top horizon- 
tal panel, scaled total kinetic energy, second horizontal 
panel, and average size of particles (d), third horizon- 
tal panel. All the measurements are done for the typi- 
cal configuration systems with 2000 particles and in zero 
gravity, but with four different boundary conditions. The 
first column of graphs, (a), shows the results for the sys- 
tem with oscillating base and glued particles on the top 
wall. The second column, (b), shows the results for the 
system with oscillating bottom wall but without glued 
particles. The last two columns show the results for the 
system without oscillating bottom wall with (c) and with- 
out (d) glued particles. When oscillations are present, 
the vertical y coordinate of the bottom wall is given by 
y = Asin(u>t), with amplitude A = (d) and frequency 
u> = u>t = 230 rad/ sec. 

In all four cases the shearing wall is rotating with 10 
rad/sec. This sets the linear shearing velocity range be- 
tween 80 (d)/sec (close to inner cylinder) and 120 (d)/sec 
(close to outer cylinder) . The effect of different shearing 
velocities is discussed in Sec. IIII CI and also later in Sec 

EH 

First we consider results without oscillations. Fig- 
ure Ufa) shows that in the case of glued particles on the 
top wall, the system reaches a state characterized by a 
significant shear throughout the domain. By comparing 
BJa) and 0Jb) , we see that there is clearly a strong ef- 
fect of glued particles on increasing shear throughout the 
domain. We emphasize that the top wall without glued 
particles, Fig.QJb), does not lead to significant shear, al- 
though the wall itself is very frictional and inelastic, see 
Tabled 

Figure 0Jc-d) shows the results with oscillations. The 
addition of oscillations accounts for a considerable slip 
velocity at the bottom wall. Similarly to the case with- 
out oscillations, the presence of glued particles greatly 
reduces the slip velocity at the top wall. Oscillations are 
also accountable for size segregation: The first two bins 
in cases (c) and (d) are populated mostly by small diam- 
eter particles, as it can be seen in the bottom panel. We 
also note (Fig. 01 top panel) a very small volume fraction 
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(*0 y bin number (b) 

FIG. 3: (a) y bin number in measurements of velocity and volume fraction profiles, (b) r bin number for wide configurations. 



close to the bottom wall. This is due to the fact that 
when the bottom wall moves downward, the adjacent 
granular layer does not follow it, except possibly small 
diameter particles. This observation is also confirmed by 
computer visualization of the granular layer. 

The distribution of volume fraction as a function of 
y is shown on the same plots as velocity distribution. 
The steady states, marked by filled circles, show that the 
volume fraction is not uniform. For each of four cases we 
observe the maximum local volume fraction « 58% for 
some intermediate y's, with dilation effects close to the 
walls, in particular the top one. 

From energy plots we estimate the equilibrating time 
for all four cases: t eq (a) — 15 sec, t eq (b) — 15 sec, 
t eq (c) = 2 sec, t eq (d) = 6 sec. These times are shorter 
for the systems with oscillating bottom wall, due to in- 
creased collision rate. Also, oscillations lead to modula- 
tion of the scaled energy; this is shown in the insets in 
the energy plots. In addition, in the systems with oscilla- 
tions, the glued particles provide an increase in collision 
rate and decrease in equilibration time; their effect on 
equilibrating time is weak in the system without oscilla- 
tions. 



B. Wide configurations with inelastic and frictional 
sidewalls 

Next, we consider the velocity and volume fraction pro- 
files for wide configuration, characterized by 5000 parti- 
cles, and 40% solid volume fraction, using typical values 
of all other parameters. The results for the same four 
boundary conditions as above are shown in Fig. [3J 

The top panel shows the velocity and volume fraction 
profiles. Qualitatively these results are similar to the re- 
sults for typical configuration shown in Fig. 0] (top panel), 
although there are some differences. For example, we 
see that the velocities are larger for wide configurations. 
This can be explained by realizing that the ratio of the 
areas of dissipative sidewalls to the area of shearing wall 
is smaller. Therefore, the influence of side walls is weaker 
in wide configuration compared to typical configuration. 
This effect also leads to larger slip velocity at the bottom 
wall. On the other hand, the volume fraction profiles are 
very similar in these two configurations. 

The second panel of Fig. [5] shows the velocity profiles 
averaged over cylindrical shells of increased radius [see 
Fig. EJb)] . To obtain these results we split the volume 
of the cell into five concentric cylindrical shells of same 
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with glued particles 
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time (sec) 




y bin number 

(a) 



without oscillations 
without glued particles 
I 



with oscillations 
with glued particles 
1 




y bin number 

(b) 



y bin number 

(c) 



with oscillations 
without glued particles 
1 
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123456789 10 
time (sec) 




y bin number 
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FIG. 4: Typical configurations results. First Panel: scaled velocity profiles (filled triangles) and volume fraction profiles (filled 
circles) in steady state regime together with the profiles in transient regimes (dashed lines). Second panel: scaled energy plots. 
Insets in energy plots: kinetic energy plotted on a shorter time scale showing the frequency of bottom wall oscillations. Third 
panel: average diameter of particles, (a), (b) - results without oscillating bottom wall; (c), (d) - with oscillations; (a), (c) - 
shearing wall with 100 glued particles; (b), (d) - without glued particles. In all simulations the dissipation parameters are as 
shown in Table 



thickness. Each shell is given a number 1 to 5 with higher 
number corresponding to the larger radius: The profiles 
for shell 1 are not shown since the reduced volume frac- 
tion near inner cylinder (see the third panel of Fig. 0) 
does not allow for good statistics. For comparison, the 
results for the velocity averaged over the whole volume 
are also replotted as solid lines. These results indicate 
that, as we change the location of the averaging shell in- 
side the cell, the changes of the velocity profiles shape 
are small, if not negligible. We note here that these re- 
sults are obtained assuming side walls are characterized 
by relatively small friction and dissipation. Therefore, 
there is no observable slow-down of the particles in the 
vicinity of the side walls. 

Figure [S] third panel, shows the radial volume fraction 
distribution obtained from averaging the volume fraction 
in each cylindrical shell. As in the case of volume frac- 
tion profiles in the vertical direction, the dashed lines 
mark the transient states and the solid lines with filled 
circles mark the steady state distributions. We observe 
significant dilation near the inner wall. This is due to cen- 
trifugal effect: Particles colliding with the top wall gain 
the momentum that, on average, has a direction along 
the tangent line to the trajectory at the contact point. 



Following this direction, particles get closer to the outer 
wall. When shearing is strong, as in the cases (a), (c), 
and (d), we also observe dilation band close to this wall. 

Finally, Fig. [SJ forth panel, shows the scaled energies 
obtained similarly to the results in Fig. ^ (third panel), 
but averaged over longer time periods to reduce noise 
(therefore, called "mean"). Comparing energy plots in 
Figs. [JJand 01 we see similar estimates of the equilibrat- 
ing times for all the boundary conditions, except the one 
without oscillations and without glued particles. Here we 
notice huge increase of equilibrating time. An explana- 
tion for this increase is that the time period to establish 
steady radial volume fraction distribution is very long, as 
illustrated in Fig[5jb), third panel. 



C. Effect of dissipation parameters and shearing 
velocity in the systems without oscillations 

Next, we discuss system's behavior when we modify 
the dissipation or the driving parameters. As a basis, 
we use the 40% volume fraction typical configuration 
system without oscillations and without glued particles, 
Fig- Bib), with one modification: to enhance the shearing 
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FIG. 5: Results for wide configurations. First Panel: scaled velocity profiles (filled triangles) and volume fraction profiles (filled 
circles) in steady state regime together with the profiles in transient regimes (dashed lines). Second panel: velocity profiles 
measured in cylindrical shells (as discussed in the text) scaled by shell's shearing velocity. Third panel: radial volume fraction 
profiles. Forth panel: scaled energy plots. As in Fig. 21 (a) and (b) are without oscillations, (c) and (d) are with oscillations, 
(a) and (c) are with glued particles on a top wall, and (b) and (d) are without glued particles. 



throughout domain, we replace the dissipative side walls 
by completely elastic side walls. 

There is a significant, qualitative difference of the ve- 
locity profiles for the system with elastic and smooth side 
walls shown in Fig. compared to inelastic and rough 
sidewalls shown in Fig.^b). Effect of shearing is much 
stronger in the systems with elastic side walls. In this 
case, the velocity profile is almost linear with very large 
slippage velocity at the bottom wall, and almost no slip- 
page at the top wall. We refer to this shape of velocity 
profile as linear -asymmetric. The asymmetry is the his- 
tory related effect: If we shear the same initial configura- 
tion by rotating the top wall in one direction and bottom 
wall in the opposite one, but with the same magnitude of 
velocity V/2, we obtain linear- symmetric velocity profile 
with equal slippage velocity at the top and bottom. 

Figure shows that the velocity and volume fraction 



profiles are very similar as shearing velocities are varied. 
Additional simulations have shown that the velocity and 
volume fraction profiles are almost identical for all the 
tested values of V between 0.1 rad/sec and 35 rad/sec. 
The energy plots (bottom panel in Fig. |SJ) show, how- 
ever, that the time to reach the steady state (equilibrat- 
ing time), t e , does depend strongly on V. This effect is 
shown in more detail in Fig. [7| which shows that t e is in- 
versely proportional to V. Therefore the same amount of 
strain is needed to reach the steady state for all explored 
shearing velocities. 

In the next set of simulations we investigate the effect 
of horizontal walls properties in the same system. Out 
of three parameters that characterize the properties of 
horizontal walls, see Table [U we vary one at a time and 
keep the others fixed. Table ITU summarizes our findings 
by describing the shearing regime for the selected ranges 
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V = 0.1 rad/sec V = 1.0 rad/sec V = 35.0 rad/sec 




1000 2000 3000 100 200 300 400 500 5 10 15 20 

time (sec) time (sec) time (sec) 



(a) (b) (c) 

FIG. 6: Linear-asymmetric velocity profiles for the case of elastic sidewalls, no glued particles, and no oscillations. Shearing 
velocity is (a) V = 0.1 rad/sec, (b) V = 1.0 rad/sec and (c) V — 35 rad/sec. 
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FIG. 7: Equilibrating times t e in simulations (open circles) are evaluated from corresponding scaled energy plots and plotted 
against shearing velocity (elastic side walls, no glued particles, and no oscillations). 



of the varied parameter. 

We find that e w and /3q w have relatively weak effect 
on the velocity profiles. However, the coefficient of fric- 
tion, fi w , influences the flow strongly. There exists a 
critical value fi^, below which the frictional force from 
the shearing wall cannot excite the granular flow. For 
our configuration this value is \i c w <~ 0.4, for the tested 
range of Vs. In what follows, we explain the influence 
of y, w on these profiles in some more detail. 

Figure |S| shows the results for fi w in the interval 
[0.4,0.9]. Here, we observe that the systems first reaches 
a metastable state of slow shearing [open symbols in 
Fig. [S] (a) and (b)]. This metastable state is charac- 
terized by negligibly small slippage velocity at the bot- 
tom wall and very large slippage at the top, as well as 
by high volume fraction close to the bottom wall. As 
time progresses, the shear of the granular particles is 
increasing very slowly. However, after certain time the 
system jumps into a stable state of fast shearing with 
asymmetric-linear profile, also shown in Fig. [S] (a) and 
(b). We call this behavior delayed dynamics. Figure |5^c) 



shows the scaled energy plots for four values of fric- 
tion which illustrate this jump from metastable to stable 
state. 

Figure ffid) shows the crossover transition time, t c . 
This transition time t c (obtained in the manner dis- 
cussed below) is defined as the time when metastable 
state ends, i.e. when the granular layer starts to pick up 
the energy fast. We find that the inverse square function 
t c = a/(p w — n%) 2 is a good representation of the depen- 
dence of t c on /j, w ; this fit is also shown in Fig.[Sfd). The 
critical friction for the onset of shearing [f w is now the 
fitting parameter of the function. Using a fitting routine, 
we obtain fj,^ = 0.33. 

Figure[5]illustrates the manner in which t c is obtained. 
This figure shows scaled energies with the tangent lines 
(shown as skewed broken lines) at the points where the 
energy increase rate is maximum. The intersection of this 
tangent line and the time axis gives an estimate of t c , and 
the intersection of the tangent line and the horizontal line 
at the level of the average scaled energy in the steady 
state gives an estimate of equilibrating time, t e . In all 
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TABLE II: Shearing regimes for different horizontal walls properties. 
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FIG. 8: Delayed dynamics regime. The shearing velocity in all cases is V — 10 rad/sec. (a) and (b) show the typical velocity 
and volume fraction profiles in metastable and stable states, using fx w — 0.45 at t — 30 sec (metastable state) and t — 150 sec 
(stable state), (c) Scaled energy as a function of time for four values of fi w . (d) Transition times, t c , versus friction: open 
circles are the simulation results, and the solid line is the best fit using inverse function t = a/(fjb — [i c ) 2 with fitting parameters 
a = 1.07 and fi% = 0.33. 
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FIG. 9: Close-up of "transition interval" allows to determine the transition time t c and equilibrating time t e , as discussed in 
the text. 



four cases, shown in Fig. EI <^e —t e — t c is approximately 
(18 ± 4) sec. This result is consistent with the result for 
equilibrating time for the same V = 10 rad/sec and for 
high friction horizontal walls (fi w = 0.9), shown in Fig.0 
For such a large fi w , t c m 0, and 5t e ~ t e . 

Next, we consider influence of shearing velocity for 
fixed friction coefficient \i w = 0.45. Figure ITUT a) shows 



the scaled energy in the transition time interval for the 
system in delayed dynamics regime for three different 
shearing velocities, V = 5 rad/sec, V — 10 rad/sec, and 
V = 20 rad/sec. This figure confirms that the transition 
times t c depend very weakly on V . They are all slightly 
scattered around the value t c 80 sec. However, the 
St e 's are following the inverse rule shown in Fig. (note 
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FIG. 10: Scaled energies in the transition interval for (a) fi w = 0.45 with three different shearing velocities, (b) fx w = 0.45 
fixed for the top wall, variable fi w for the bottom wall, fixed shearing velocity V = 10 r ad/ sec. 



that t c — for fx w = 0.9 used in Fig.Q). Therefore, <5t e 's 
strongly depend on the shearing velocity, while being al- 
most independent of the coefficient of friction. 

To understand precisely the role which the horizon- 
tal walls play in the mechanism of delayed dynamics, 
we have performed simulations with fixed coefficient of 
friction of the top wall n w (top) = 0.45 and varied coeffi- 
cient of friction of the bottom wall. The results, shown 
in Fig. HUf b). clearly indicate that fi w (bottom) does not 
have the determining effect on the transition time t c . For 
example, changing it from 0.9 to 0.25 leads to less then 
5% decrease of t c . However, the higher slope of energy 
increase for smaller friction coefficients leads to smaller 
St e . We also note that for these low /j, w (bottom) the ener- 
gies are rising to the higher values. This is an indication 
of higher overall velocities of the particles in the steady 
state due to higher slippage at the bottom wall for lower 
/j, w (bottom) . 

Based on the results shown in Figs. 171 101 the mech- 
anism of the delayed dynamics can be explained as fol- 
lows. The friction of the top wall controls the amount of 
the tangential momentum given to the particles colliding 
with it. Thus the energy input to the granular system 
is determined by the top wall friction, as well as by the 
collision rate of the particles with this wall. When the 
initial configuration is subjected to shear, the particles 
immediately dilate close to the top wall and form the 
cluster of high volume fraction close to the bottom wall. 
Reduced number of particles close to the top wall then 
significantly decreases the input of energy. If the friction 
of the top wall is low enough, below some critical value 
fj,^, all the energy is quickly dissipated without induc- 
ing shear in the system. For higher frictions, the energy 
transfered to the granular system starts to accumulate 
slowly: The shearing throughout the sample increases, 
until, at some critical time t c , the shearing between the 
dense cluster and the bottom wall is big enough to dilate 
a region next to this wall and to push the cluster closer 
to the top wall. This process is sped up if (i w (bottom) 
is smaller. As a result, the collision rate with the top 
wall increases, and the energy from this wall is rapidly 



transfered into the kinetic energy of the whole cluster, 
until the stable state characterized by linear asymmetric 
profile forms. Thus, the time span of metastable state 
t c is determined mostly by the friction of the top wall, 
and the time 5t e is determined mostly by the shearing 
velocity. 



D. Effects of oscillations 

To test the effect of the amplitude A and frequency / 
of oscillations, it is more convenient to go back to typi- 
cal dissipative sidewalks used in Sec. lIII Al The reason is 
that in the case of smooth side walls the slippage velocity 
at the bottom wall is approaching the value of shearing 
velocity (this is shown later in this Section), making the 
analysis of the velocity profiles difficult. So, as a base 
system we use the "oscillating bottom wall, no glued par- 
ticles" configuration shown in Fig. Efd). 

Figure 1111 shows the results for velocity, volume frac- 
tion and equilibrating times for the systems where: (a, b) 
/ is fixed and A is varied, (c, d) A is fixed and / is varied. 
The characteristic feature of the systems shown in Fig.llll 
is the high volume fraction (v rs 50 — 60 %) band of the 
granular layer in the middle of the cell height, confined 
between two thin layers of low v close to the bottom and 
top walls. Inside this band the local shear rate is rather 
low, compared to the overall shear rate V/H. The high- 
est local shear rates are observed only close to the top 
wall, while the high slippage velocities are observed at 
both top and bottom wall. 

Figure lllf a) shows the amplitude dependence of the 
velocity and volume fraction profiles for fixed frequency 
/ = 36.6 Hz. We observe that the slippage of the bot- 
tom wall depends strongly on A, with higher A's leading 
to more slippage. Also, as A increases, the peak in the 
volume fraction profile moves further away from the oscil- 
lating wall. The equilibrating times, Fig. HiT b'). decrease 
with an increase of A in a roughly linear manner. 

Figure lllf c) shows the frequency dependence of the 
velocity and volume fraction profiles. Here, A is fixed 
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FIG. 11: Velocity and volume fraction profiles for (a) fixed frequency, ui = uj t (/ = 36.6 Hz), and different amplitudes of 
oscillations; (c) fixed amplitude, A = (d), and different frequencies. Also, equilibrating times are shown as a function of (b) 
amplitude and (d) frequency of oscillations. The volume fraction symbols follow the line pattern of the amplitude (a) or 
frequency (c) ones. 



to 1 (d). The results show, again, that an increase of 
the intensity of oscillations, i.e. frequency, leads to an 
increase of velocity throughout the system. However, the 
volume fraction profiles change very little with /. The 
equilibrating times, Fig- IHf d). decrease with an increase 
of /, with the lower rate of decrease for higher /'s. 

Previously, by comparing Fig. 0Jb) and Fig.H3 we saw 
significant effect of sidewall properties on velocity and 
density profiles for the systems without oscillations. Fig- 
ure 1121 shows the corresponding results as properties of 
the sidewalk are modified in the case of oscillating bot- 
tom wall, and absence of glued particles. To simplify the 
comparison, Fig. ll2f a) shows again the result for the typ- 
ical case, Fig. Eld). Setting ^ s = 0, e s = 1 [Fig. \T%b)} 
makes almost all particles move with shearing velocity. 
Setting fi s — 0.5, e s — 0.6 [Fig. Et c )] creates already 
enough dissipation of energy to almost completely re- 
sist the shearing force: even oscillations cannot improve 
shearing when sidewalls are very frictional. Next we ex- 
plore whether it is /i s or e s that influence the velocity 
profiles. Therefore, in Fig. I12f d^ we show profiles and 
scaled energy data for [i s = 0.1, e s = 0.6. By comparing 
(a) and (d) we see that the elasticity of the side walls 
does not have the observable effect on the velocity pro- 
files. Therefore, fj, s plays the crucial role, similarly as 
observed regarding \i w in Sec. lIII CI The volume fraction 
profiles and equilibrating times, on the other hand, very 
weakly depend on the side wall properties. 

IV. NORMAL AND SHEAR STRESSES ON 
THE BOUNDARIES 

Stresses in granular systems have been the focus of 
many studies, because of their importance in a num- 
ber of engineering designs. A significant part of these 
studies22i22i2£ deals with theoretical continuum models 
where stress is considered as a mean local quantity. 
Therefore, the validity of these models depends on the 
strength of stress fluctuations on the scale which defines 



"locality" . The knowledge of stress distributions becomes 
crucial here. There are many experimental, numerical, 
and theoretical studies of stress fluctuations ^2&2L22i22i 
however, these studies deal with high volume fraction 
of static or slowly sheared granular systems. In these 
dense systems the distribution of stresses is character- 
ized by an exponential decay for large stresses. It is 
speculated^ that the exponential tails in stress distri- 
butions are related to the presence of force chains in 
these systems. However, recently, Longhi, Easwar and 
Menon^S reported the experimental evidence of exponen- 
tial tails in stress distributions in rapidly flowing granular 
medium, i.e. in the system where it was unlikely to find 
force chains in the traditional "static" sense. 

The results presented in this Section deal with the 
stress distributions at the physical boundaries of rapidly 
sheared granular systems that are not dense, i.e. 40% 
volume fraction. Before presenting these results we 
should note that in order to make maximum connection 
with existing experimental results and theories, we calcu- 
late the stresses in a manner that is very similar to experi- 
mental methods, i.e. we also introduce in our simulations 
stress sensors on the boundaries (more details provided 
below). In the low volume fraction and/or rapid flow 
regime these sensors can report zero stress during some 
measuring intervals, see, for example, Ref.l30l Therefore, 
the total distribution of stresses must contain the distri- 
bution of non-zero stresses and a weighted delta function 
to account for zero stresses. In what follows we make a 
clear distinction between zero and non-zero stresses and 
extend all relevant theoretical arguments to account for 
the presence of zero stresses. 

We present the results for the stresses on the bound- 
aries for the typical configuration systems shown in 
Fig.^J More detailed analysis will be given for the partic- 
ular system without oscillations and with glued particles, 
Fig. 0fa). Of particular interest here are the distribu- 
tions of normal and shear stresses on the top and bottom 
walls. To calculate these stresses we cover the respective 
boundaries with a grid of square stress sensors, as shown 
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FIG. 12: Effect of sidewall properties on velocity and volume fraction profiles, (a) Base system, e s 
(b) Completely elastic and smooth, (c) e s = 0.6, fj> a = 0.5, flo s = 0.35. (d) e s = 0.6, fi s = 0.1, /3q s - 
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FIG. 13: Stress sensor areas between inner and outer cylinder of Couette cell, (a) Square sensors. Only the data from the 
sensors completely inside the volume of the cell are used, (b) Sector-shaped sensors. 



in Fig. I13f a) . All sensors have the same area A s = d 2 s . 
Typically, we set d s — 1.09 (d), A s = 1.19 (d) 2 , and ana- 
lyze the data only from the sensors that are completely 
inside the volume of the cell. Later in this Section we use 
sector shaped sensors, Fig. I13f b). which are convenient 
to set large sensor areas. The instantaneous normal and 
tangential stresses on a sensor are defined in the following 
way: 



CTt 



Here i counts all the collisions between the particles and 
the walls in the sensor area during a time interval At. We 
refer to this time interval as averaging time. Ap n (i) and 
Apt(i) are the change of particle's normal and tangential 



components of momentum, p s is the density of the solid 
material, di is the diameter of a particle participating in 
the collision i, and vectors v and v' are the velocities 
of the colliding particle before and after a collision, re- 
spectively. Subscripts n and t refer to the normal and 
tangential components. Our typical choice for the time 
interval is At = T w /10. In order to get sufficient data 
for the stress Probability Distribution Function (PDF), 
we measure instantaneous stresses on all sensors every 
At for total of at least AT = 1000 At. For the con- 
sidered geometry, Fig. I13f a). we have the total of 124 
sensors between the inner and outer wall. When present- 
ing the results for stress distributions we show the scaled 
histograms (PDFs) of locally measured normal or tan- 
gential stresses, a n and at, as a function of a n /(a n ) or 
a t /{a t ), respectively. The averages are defined as 



T (n,t), 



= (J2&P(n,t)(i))/A base /AT 



(10) 
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A base = ir(R 2 - R 2 ). (11) 

Here i counts all collisions between the particles and the 
top or bottom wall during the whole measuring period 
AT = 1000 At. 



A. Normal stress distributions 

Figure 1141 shows the results for normal non-zero 
stresses on the top and bottom walls in the typical config- 
uration for the same four boundary conditions that were 
discussed in the preceding Section (see Fig.^J. The main 
features of the stress distributions are the increase from 
zero for small stresses and the decrease, in most cases ap- 
proximately exponential, for large stresses. We note that 
it is well established experimentally that the exponential 
tails are the characteristic features of the stress PDFs in 
static granular materials at volume fractions close to the 
random close packing. ?L2£i2I Although the stresses we 
measure are not static, nevertheless, the approximately 
exponential tails are observed in all cases shown in Fig. 1141 
(except perhaps for very large stresses). This is evident 
from the best fit of the following functional form 

PnM =ca"e^ /CTo ; a = a /{a), (12) 

shown in Fig. El as solid lines; we use (a) = (a n ) as 
defined in (|10|l . The subscript nz is used to emphasize 
that the distributions are build from non-zero stresses. 
This is important because certain fraction of sensors reg- 
isters zero stress due to the fact that the volume fraction 
(typically 40%) is not sufficiently large to warranty that 
there is always a particle in contact with a sensor. In- 
cluding zero stresses in the definition of P would prevent 
us from formulating the fitting function (|12fl and relate 
our results directly to experimental and theoretical re- 
sults. We discuss the consequence of this approach in 
more details below. 

Form (|12f) can be thought of as the generalization of the 
theoretical predictions for static stress distributions 42 . 
For example, the original q model proposed by Liu, 
Nagel, Schecter, Coppersmith, Majumdar, Narayan, and 
WittenSi^ and extended later by other authors f^ 3 ^ pre- 
dicts the distribution (|12fl with n = N c — 1, where N c is 
the number of force transmitting contacts between a par- 
ticle in one layer and particles in the adjacent layer. In 
three dimensions and fee close packing N c — 3, so n = 2. 
Using different approach, Edwards and Grinev 35 predict 
the distribution (|12fl with n = 1/2. n can also depend 
on the number of contacts between sensor and particles 
during averaging time, as will be shown in Sec. IIV A 21 
Therefore, it makes sense to think of n as a fitting pa- 
rameter. Other parameters, c and o~q, can be expressed 
in terms of n and the average stress (a) nz (measured di- 
rectly in simulations) using the following expressions: 

/>oo 

1 = / P nz (a)da = c[T(n + l)^ 1+n) ] (13) 
Jo 



{a)nz = \ aP nz (a)da = a (n+l) (14) 
Jo 

In Fig. E] we show the fit of the functional form 1121) 
with n as the only fitting parameter. Other parame- 
ters, calculated using l(T3"|) and (fHIJl . and average quanti- 
ties, obtained directly from simulations, are also shown. 
These results show that the average stresses (a) (indi- 
cated in the horizontal axes labels on each plot) at the 
top and bottom walls are equal within the statistical er- 
ror. However, the peak value and the width of a dis- 
tribution is different in each case. In particular, for the 
systems without oscillations the distributions are wider 
at the top wall, while if oscillations are present, they are 
wider at the bottom wall. Before discussing this result, 
we first consider the effects of volume fraction, collision 
rate, averaging, and correlations on the properties of the 
stress distributions. 



1. The width of stress distributions 

In previous section we have shown that a distribution is 
determined by the parameter n and the average stresses 
{o~)nz and (a). Here we define the width of a distribution 
in terms of the above quantities and discuss its depen- 
dence on the various properties of granular system. 

Let the width of a distribution be defined as the root 
mean square of all measured stresses (including zero 
stresses) divided by the mean stress 

width = vwtew _ vwtew _ y/m ^ l 

[ (T ) (°7 

(15) 

(here we use (a) = 1). The distribution itself takes the 
form of a linear combination of the distribution of non- 
zero stresses (|12l) and the delta function to account for 
zero stresses 

P(a)=C 1 P nz (d-) + (l-C 1 )S(0) (16) 

The constant C\ signifies the fraction of non-zero stresses 
in the total distribution. Then, the mean and the mean 
square of the stress are as follows 

/>oo 

{a) = / aP{a)da = C 1 {a) nz = 1 (17) 
Jo 

/•OO 

(a 2 ) = / a 2 P{?r)da = C l0 -l{n + \){n + 2) (18) 
Jo 

From the first equation we obtain C\ = l/{cr) nz and, 
using (|14l) . we arrive at the following expression for the 
width of distribution: 

width = ^/a {n + 2) - 1 = v^o + WU* ~ 1 (!9) 

Expression l|19|) reduces to the simple ^Joq — (n + 1) -1 / 2 
in the case of absence of zero stresses when all sensor reg- 
ister at least one strike per At and (a) nz = (a) = 1. In 
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FIG. 14: Log plot of PDF for normal non-zero stresses on the top wall (top panel) and on the bottom wall (bottom panel) in 
typical configuration with four boundary conditions: (a) and (b) without oscillating bottom wall, (c) and (d) with oscillating 
bottom wall, (a) and (c) with glued particles, (b) and (d) without glued particles. Average normal stress per area (<r) is in 
units of p s {d} 2 /sec 2 . Solid lines: the best fits of the function 111 2t . with the values of the fitting parameter n and corresponding 
cro and (a)nz shown on each plot. 



this regime, the width is determined by averaging effects 
and correlations, as we discuss later. On the other hand, 
in the case zero stresses are present, we have a significant 
dependence of the width of a distribution on the average 
non-zero stress (a) nz . In a limiting case (a) nz 3> 1, ap- 
propriate to most systems shown in Fig.^J equation l|19|) 
reads: 



width 




1 « V(°)n Z (20) 



where we have used the fact that for (a) nz 3> 1 sensors 
do not register more then one particle at a time and for 
single particle contacts n — constant, as explained in next 
Section. Equation (|20|l is consistent with visual estimate 
of the widths of the distributions in Fig. [21 The larger 
(ff)„ z (shown in the right top area on each plot), the 
wider is the distribution. Next we show that {d) nz in 
turn depends on the collision rate, area of the sensors 
and averaging time. 

Let N t = n s riT be the total number of stress mea- 
surements, obtained from all n s sensors during «t time- 
intervals. Out of N t stresses, N nz < N t are non-zero. 
The non-zero average stress depends on the size of sen- 
sor A s and averaging time At, 



N t = 1 
N nz w A s At 



(21) 



where w is the frequency of non-zero stress events per 
unit area. In the case of low collision rate with the sen- 
sors, when no more then one particle strikes A s during 



At, w has a meaning of particle-wall collision rate per 
unit area. 

Relations (|20|l and l|21|) allow us to explain the distri- 
butions in Fig. ^31 The sensor areas and averaging times 
are fixed so it is the rate of non-zero events w which de- 
termines the average and width of a distribution. To the 
first approximation w oc the volume fraction close to the 
respective boundary. From Fig. 01 (top panel) , in the sys- 
tems without oscillations the volume fraction is smaller 
at the top wall compared to the bottom wall. Therefore, 
the distributions are wider at the top compared to the 
bottom. In the systems with oscillations the volume frac- 
tion is smaller at the bottom wall compared to the top 
wall, therefore, the distributions are wider at the bottom. 



2. Averaging effects 

In this section we study one of the factors that can 
affect the parameter n, i.e. the number of contacts, N, 
between the particles and a sensor during the averaging 
time At. 

Let us define the parameter uq that is characteristic for 
the distribution of stresses due to single particle contacts: 



N = 1 



(22) 



If more then one particle collide with the sensor of area 
A s during time At (N > 1), then the stress distribution 
is different due to averaging effects. For instance, when 
two particles strike the same sensor during At, N = 2, 
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the probability of registering total stress F depends on 
the probability of one of them contributing the stress 
a < F and another one the stress F — a. Assuming both 
considered particles have independent distributions (|22p. 
the distribution of their collective stress F is Pi(F} — 

J P(a)P(F — a)da. Using and the integral identity 

f F 

/ a r (F-a) q da = F( r+q+ V r(r+l)T{q+l)/r{r+q+2) 
Jo 

(23) 

valid for any positive real numbers r and q, and in par- 
ticular for r = q = riQ, we arrive at the following distri- 
bution for the stresses generated by double strikes: 

P 2 (F) = c 2 F ( - 2no+ ^e- p ^«- 1 F = F/(a) (24) 

where c 2 is a normalization coefficient. Applying the 
same argument to the case when N particles strike the 
sensor during At we obtain the following distribution of 
stresses: 

P N {a) = c 5-W"o+ 1 )- 1 ) e -*/° n (25) 

The above calculation sets the relation between the pa- 
rameter n and the number N as n — N(no + 1) — 1, 
see lfT2"|l . Therefore, under the assumption that parti- 
cles are not correlated we obtain that the distributions 
of stresses generated by multiples collisions are charac- 
terized by higher power n. We will see below that this 
assumption is justified for the results shown in Fig. 1141 

This conclusion may be used to explain the increased 
values of n observed in the distribution of stresses on 
the bottom wall in the systems shown in Fig. I14f c') and 
Fig. I14f d). In these systems the bottom wall is vibrated 
and most of the stress data are collected during the phase 
when the bottom is rising up. This phase is definitely 
characterized by increased compaction of particles close 
to the sensors and an increased number N. 

However, an increase of N cannot be responsible for 
the large n's observed on the top wall in the system shown 
in Fig. H4T b). and, to lesser extend, Fig. H"4T a). In par- 
ticular regarding Fig. ll4f b). we expect that the source of 
large n lies in the separately verified result that the typ- 
ical normal components of the velocities of the particles 
colliding with the top wall is large, leading to different 
stress distribution. We note, however, that the values of 
n given in Fig. I14f a-b) are not as accurate as the rest of 
n's shown in Fig. 1141 due to small volume fraction there 
(see Fig. |3Ja-b)). Therefore, the number of collisions is 
small - e.g., we count approximately 500 collisions dur- 
ing the whole time span of the simulations presented in 

Fig. mb). 

The results of Fig. El allow us to estimate the param- 
eter uq defined earlier in this Section as the value of n 
corresponding to the distribution of stresses registered by 
the sensors with N = 1. Because uq < n (see J250). the 
smallest found n provides an upper bound on n Q . The 
lowest values of n are found in the distributions shown in 
Fig- El ( a ) bottom, (b) bottom, (c) top and (d) top. For 



these cases n < 1, therefore no < 1. This result is very 
different from no — 2 in the g-model: this should be no 
surprise because g-model assumes high volume fraction 
and static configuration. 



3. Correlations 

The dependence of the normal stress distributions on 
the number of contacts per sensor can be described by 
(|25|l only when the particles participating in contact with 
a sensor are not correlated. However, there are cer- 
tain regimes and conditions of granular flow when this 
assumption may not be correct. For example, Miller, 
O'Hern and Behringer— studied the response of the force 
distribution to the number N, controlled in their case by 
the size of particles, and found that in dense granular flow 
the distribution widths were approximately independent 
of N. This independence was attributed to the presence 
of correlations due to presence of force chains. 

In what follows we present the systematic study of 
the effect of sensor size (and, hence, the number N) on 
the bottom wall stress distributions for our typical con- 
figuration without oscillations and with glued particles, 
Fig.QJa). This study will allow us to confirm the validity 
of our predictions (|19[) and (|25|) in a more quantitative 
way and to check for possible correlations. 

In order to explore large range of sensor sizes, here 
we use sector-shaped sensors, obtained by dividing the 
bottom wall area by a number of radial lines drawn 
from the center of the cell toward outside boundary, see 
Fig. I13f b). Each sensor has the same radial dimension 
equal to the distance between the inner and outer cylin- 
der, and the same angular dimension, set to the values 
2ir/m, m G [2, 252]. Therefore, the area of sector-shaped 
sensors varies between 1.0 and 125.6 (in units of (d) 2 ), 
with the largest sensor being half of the bottom wall. For 
selected values of sensor areas we have confirmed that the 
results do not depend on the sensor shape. We note that 
all the results presented so far were obtained using small 
sensors of the area 1.19 (d) 2 . 

Figure I15f al shows the results for width of distribu- 
tion (|15|) as a function of A s for fixed values of 5t in the 
range between 7^/100 = 0.0002 sec and 20T lu = 0.5463. 
This range includes our typical value At — T w /10 = 
0.0027 (open triangles) used to obtain the distributions 
in Fig El To help interpretation of the results, we plot 
in Fig. EJb) the value of (cr) nz , in Fig. I15f c) the aver- 
age number of contacts per sensor, and in Fig. I15f d~) the 
widths of distributions versus the product of the averag- 
ing time and the sensor area. 

This figure shows that if zero stresses are present, i.e. 
(&)nz > 1) the widths scale with the sensor area as A~ 05 , 
which is the scaling predicted by This result applies 
to all sensor sizes when At = 0.0002 sec, to the sensors 
of A s < 20 when At — 0.0027 sec, and to the sensors of 
A s < 2 when At = 0.0273 sec, see Fig.lTsTb). 

When zero stresses are not present ((a) nz = 1) the 
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FIG. 15: (a) The widths of force distributions as defined by l)15[l versus the sensor area A a (in units of (d) 2 ) for six values of At 
(shown below the plot), (b) Average non-zero stress {a)„ z as a function of A s for the same six averaging times (for At > 0.2731 
(d')nz = 1). (d) Average number of contacts per sensor as a function of A 3 . (c) Collapse of the width data against the product 
At A„. The solid lines have specified slopes. 



width Ijl9|l reduces to: 

width = V^o = (N(n + l))-°- 5 (26) 

Here N oc A s is a number of contacts per sensor from H25|) 
and absence of correlations is assumed. This relation l|26|l 
is consistent with the Aj - 5 in Fig. ll5f a) in the case At = 
0.0027 sec and A s > 20. However, when At > 0.0273 (for 
these large At's, the zero stresses are absent for almost 
all considered sensors, see Fig. HHT b)) the widths decrease 
with an increase of A s at a slower rate then predicted by 
(|2*tj|l . We observe the slowest decrease, approximately 
Aj ' 21 , for the sensor sizes 10 < A s < 100, and for the 
averaging times of 0.2731 sec. 

We expect that the origin of this slower decrease of 
the widths is the presence of correlations between the 
particles striking the same sensor (which may include 
self-correlation, i.e., the same particle colliding with a 
same sensor on multiple occasions during a given At). 
For example, the distribution of stresses due to N com- 
pletely correlated particles acting on each sensor would 
have n = uq, leading to the widths independent of A s . 
Therefore, the stronger are the correlations between par- 
ticles, the lower is the exponent in the width versus A a 



function. In Fig. I15f a) we can see that stronger correla- 
tion effects (i.e. slower decrease of the widths with an 
increase of A s ) occur for longer At and larger A s . This is 
even more evident when the width is plotted against the 
product At times A s , see Fig. I15f d): the slope changes 
gradually from —0.5 for small At A s to about —0.3 for 
large At A s . The collapse of almost all data points on ap- 
proximately single curve in Fig. I15f d) indicates the fact 
that the same averaging or correlation effects result if ei- 
ther At or A s is increased by the same factor. Possible 
exception to this rule are very long At's, which we discuss 
below. 

We note that Miller, O'Hern and Behringer, 26 have 
found evidences of strong correlations when they es- 
timated experimentally the width cx A® for At = 
0.0005 sec and 5 < A s < 80. However the overall vol- 
ume fraction of their samples was considerably higher 
(because of the heavy top wall supported in gravitational 
field by the granular material), therefore the spatial cor- 
relations can be explained if one assumes the normal 
stress on the bottom wall is applied through the network 
of force chains*^ The results presented here are obtained 
for smaller volume fraction and in the absence of gravity. 
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However, in our case we expect that spatial correlations 
may arise from statistical fluctuations of the local volume 
fraction. These fluctuations lead to occasional formation 
of the denser structures of particles, stretched from the 
top to bottom wall: the presence of friction between the 
particles insures the lifetime of these formations is rel- 
atively long. As in the case of force chains, the dense 
structures carry more stress then the "free" particles in 
the "inter-structure" space, thus introducing spatial cor- 
relations. Unlike the force chains, the dense structures 
do not form dense force network and are not present al- 
ways. They may form and disappear allowing for longer 
"structure- free" periods. This concept is consistent with 
the fact that our correlations are weak and the correla- 
tion regimes are detected only when At and A s are large 
enough to include such a structure event. In addition the 
average number of contacts with sensors, see Fig. I15f c'l , 
is larger then 10 in the regimes where (|26H fails. This 
fact suggest the possibility of frequent returning contacts 
of the same particle with a sensor, which in event-driven 
simulations is a signature of a particle pressed to the sen- 
sor from above, as one would expect for a particle that 
is member of a dense structure. The existence of similar 
structures (named as "transient force chains") was also 
suggested in Ref. 

We note that the dependence of the width on A s be- 
comes stronger again in the case of very long At's, see 
the results for At — 0.5463 in Fig. I15f a1. Simple esti- 
mate shows that this change in trend occurs when At is 
so long that a typical particle, moving with a typical ve- 
locity, travels across a sensor in the time that is shorter 
than At, therefore decreasing this self-correlating effect. 
We currently investigate the scenarios proposed here in 
more detail and plan to present the results elsewhere. 



B. Tangential stress distributions 

Tangential (shear) stress is defined by @. Figure iPol 
shows the tangential stress distributions (filled triangles) 
scaled by mean tangential stress, a = at /(at) for the 
same boundary conditions as discussed in Fig. El For 
the convenience of comparison, the normal stress distri- 
butions from Fig. E| are also shown (open circles). 

The average tangential stress, shown in the horizontal 
label on each plot, is positive on the top wall and neg- 
ative on the bottom. This reflects the fact that the top 
wall is mostly accelerating the particles while the bot- 
tom wall is slowing them down. The absolute value of 
the average tangential stress (at) tells us the magnitude 
of torque and power needed to keep the top wall rotat- 
ing and the bottom wall stationary. It is considerably 
larger for the systems with glued particles compared to 
the ones without glued particles, since the glued parti- 
cles enhance momentum exchange between the wall and 
the free particles. We also note that \(a t )\ is larger on 
the top wall then on the bottom by a factor of 3 for 
the systems without oscillations, and by a factor of 5 



for the systems with oscillations. This factor is larger in 
the latter system since oscillations increase the vertical 
velocity fluctuations and the dissipation of the momen- 
tum on the sidewalls is increased. Since the momentum 
conservation requires that the momentum applied to the 
system through the top wall be equal to the momentum 
dissipated on both sidewalls and on the bottom wall, less 
momentum is left for the bottom wall in the systems with 
oscillations. 

Similarly to the normal stresses, we observe the ex- 
ponential tails for a > a$ where ao is the location of 
the peak of a distribution. For most considered cases, 
these tails have the same slope and location as the tails 
of the normal stresses. This suggests strong correlation 
between the normal and tangential stresses. However, for 
a < ao the distribution of tangential stresses differs qual- 
itatively from the distribution of normal stresses. The 
PDF's of tangential stresses decrease to zero in an ex- 
ponential manner, allowing for stresses in the direction 
opposite (at). However, in most cases this different be- 
havior has only weak influence on the width of the dis- 
tributions. 

An exception is the case of oscillating bottom wall 
(Fig. I16f c) and (d) bottom panel) where the tangen- 
tial distributions are noticeably wider then the normal 
stress distribution. This can be explained by noticing 
that when bottom is moving up, the particles next to the 
bottom wall are constrained by increased volume frac- 
tion to move predominantly in the horizontal directions, 
therefore increasing fluctuations of tangential stresses. 

To better understand the relation between the normal 
and tangential stresses we study in more details the time 
series of these stresses. Figure IT7I top two panels, shows 
the values of a n /(a n ) and at/ (at) as a function of time 
registered by a single sensor on the bottom wall for the 
four boundary conditions as in Fig. El Simple visual in- 
spection of Fig. 1171 top two panels, shows that the normal 
and tangential stresses are correlated. This observation 
is confirmed by comparing the normal stress autocorrela- 
tion function with the appropriate cross-correlation func- 
tion, as we discuss next. 

Consider the time signals a n (t) and at(t) shown in 
Fig. El nrs t and second panels. If these signals are cor- 
related then we can write 

a t [t) = Ka n {t) + R{t), (27) 

where R(t) is a noise term such that (R) = and 
K = (at) / (a n ) = 1. Then the cross-correlation func- 
tion defined as (a n (to)a t (to + t)) must be equal to the 
autocorrelation function defined as (a n (to)a n (to + t)). In 
the definition of these time correlation functions, averag- 
ing is performed over all initial times to. To reduce noise 
we also average the results of all sensors. 

FigureEl third panel, shows that to within small error 
(o'n.o't) — (®n&n), therefore, confirming the relation 127|) 
with (R) = 0. This figure also shows that the correla- 
tion functions decay very fast within an interval shorter 
than averaging time. For longer times, the correlation 
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FIG. 16: Log plot of PDF for tangential non-zero stresses (filled triangles) on the top wall (top panel) and on the bottom wall 
(bottom panel) in the typical configuration with four boundary conditions: (a) and (b) without oscillating bottom wall, (c) 
and (d) with oscillating bottom wall, (a) and (c) with glued particles, (b) and (d) without glued particles. Average tangential 
stress per area (er) — (a t ) is in units of p s {d} 2 /sec 2 . For comparison, normal stress distributions from Fig. ll4l are replotted here 
(open circles). 



function decreases until it reaches its uncorrelated value 
equal to unity in the systems without oscillations, (a) and 
(b), and it oscillates around this uncorrelated value with 
the frequency of driving oscillations for (c) and (d) . We 
note that larger particle velocities lead to faster decay of 
correlation functions in Fig. ll7f a) compared to ll7f b') (see 
also Fig. EJa-b), top panel). The uncorrelated value for 
both correlation functions is equal to unity because an 
average of product of two uncorrelated stresses is equal 
to the product of average stresses. 

The results of this section confirm strong correlation 
between the normal and tangential stresses. This corre- 
lation renders that many conclusions obtained for normal 
stresses also apply to the tangential stresses. 



V. CONSTANT VOLUME VERSUS CONSTANT 
PRESSURE SIMULATIONS 

All the results we have discussed so far are obtained for 
the systems characterized either by fixed constant vol- 
ume, as in the simulations without oscillations, or by 
prescribed volume dependence on time, as in the simu- 
lations with oscillating bottom wall. We refer to these 
cases as controlled volume simulations. A special case 
of controlled volume simulations are the constant volume 
( CV) simulations where the volume of the system is kept 
strictly constant (no oscillations). 

It is, however, well known, that the volume constraint 
can significantly affect the granular flow at high volume 
fractions. The system can be easily locked in a jammed 



state. That is why dense sheared granular flow experi- 
ments are commonly carried out using adjustable volume, 
but controlled pressure. Miller, O'Hern and Behringer, 26 
for example, report the results for the sheared granular 
flow in a Couette cell, where the weight of the top wall 
provides constant pressure. Also, all free surface granular 
flows in gravitational field, for example, the flow of the 
granular matter down an inclined plane fi£L are the ex- 
amples of controlled pressure scenarios. In these cases the 
weight of the granular matter itself controls the pressure 
both inside the system and on the boundaries. Recently, 
new studies emerged that indicated the importance of 
the difference between controlled volume and controlled 
stress boundary conditions. For example, Aharonov and 
Spar k q 2 i 38 use the 2D numerical simulations to study the 
response of the dense sheared granular system to the ap- 
plied boundary pressure. They find different shearing 
modes determined by pressure; These modes can have 
very similar volume fractions, but different microstruc- 
tural organization. 

In this section we describe the simulations of sheared 
granular flow in a system where the pressure on one of the 
walls is controlled at all times. We refer to these bound- 
ary conditions as controlled pressure simulations. When 
pressure is kept constant we have the case of constant 
pressure ( CP) simulations. In laboratory experiments the 
CP condition can be realized not only by setting the con- 
stant weight load on a movable top wall of a Couette 
cell 2 ^, but also by setting the movable wall, for exam- 
ple, a bottom wall, on a compressed spring, see Fig. ^| 
and Ref. Any increase of the pressure inside the cell 
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FIG. 17: First panel: normal stress signal on the bottom wall during randomly chosen interval of 2 seconds. Second panel: 
tangential stress signal on the bottom wall for the same 2 seconds interval. Third panel: stress auto- and cross-correlation 
functions. The inset in (b) shows the correlation functions on the longer time scale, (a) and (b) - without oscillating bottom 
wall, (c) and (d) - with oscillating bottom wall, (a) and (c) - with glued particles, (b) and (d) - without glued particles. 



bottom wall 




compressed spring 



FIG. 18: Schematic representation of Couette cell in vertical 
cross-section with movable bottom wall: possible realization 
of constant pressure experiment. 



results in an increase of the volume of the cell and vice 
versa. The compressed spring realization is especially 
relevant if the experiment is conducted in zero gravity. 
In numerical simulations we adopt the following algo- 



rithm (which is an approximation of a compressed string) 
to control the pressure. The whole timespan of simula- 
tion is divided into small time intervals At. At the be- 
ginning of each time interval we calculate the total stress 
on the bottom wall, er, using the impacts of the particles 
during the previous time interval. More precisely, we use 
Eq. |JSJ with A s = A base and put At = At. We then 
compute the vertical velocity of the bottom wall, Vb, as 
follows: 



Vb 



—Vb, a > a + 5a/2 

0, a - 5a/2 < a < a + 5a/2 

+Vb, a < er — 5a/ 2 



(28) 



where a and 5a are the prescribed stress and stress tol- 
erance, and Vb is a prescribed maximum instantaneous 
velocity of the bottom wall. Within each time interval 
At we run the controlled volume simulations with con- 
stant velocity of the bottom wall given by Eq. (|2*5|l . The 
resulting changes of the volume of the cell adjust the 
stress back to the prescribed value. This model leads 
to vibration of the bottom wall in steady state; similar 
vibrations are also observed in experiments^ 
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We consider the system without oscillations and with 
glued particles on the top wall, Fig. Ufa), as a basis for 
all our CP simulations. We set the constant stress on 
the bottom wall being equal to the average stress (cr n ) 
on this wall which we measured in CV simulations with 
V = 10 r ad I sec. The choice of the parameters used in 
the model J2EI is 



'• 2 

seer 

seer 
At = T w /25 

v b = (d) /At/1000 = 0.915(d) /sec 



a = 222 p, 



5(7 = 20 p ; 



(29) 



(30) 



Figure I19f a) shows the PDF for the stresses on the 
complete bottom wall in the case of CP simulations (cir- 
cles) using parameters given by l(2*9")) - lj3T)| and At = At. 
For comparison, the results of CV simulations for the 
same A s and At are shown on the same plot (triangles). 
We see that the CP distribution is wider compared to 
CV case. To understand this difference it is useful to dis- 
tinguish three sets of collected stress data, each set con- 
taining the data from one of three phases of the bottom 
wall motion: moving up, staying put, and moving down 
(this separation of the stresses is possible only because 
At = At, i.e. during the averaging time the bottom 
wall does not change its velocity). In each phase we have 
controlled volume evolution of the system with a differ- 
ent mean stress: highest in the phase of the bottom wall 
moving up and lowest in the phase of the bottom wall 
moving down. The resulting CP distribution must corre- 
spond to the superposition of the distributions from three 
phases. The degree of widening of the stress distribution 
can be estimated using the following argument. 

If Avy CV ^ is a change in the normal component of ve- 
locity of a particle interacting with a stationary sensor 

(CP) 

and Avy is a change in the normal component of the 
velocity of the same particle interacting with a sensor 



that has the velocity Vb, then Aw 



(CP) 



Av 



(CV) 



(note that Vb can take three different values, see (|28H 'I. If 
we have N particles of average diameter (d) interacting 
with a sensor A s during At then, using the definition JS}, 



r (CP) 



A* At 



(31) 



where a^ P ^ is a stress on a moving sensor, a^ V ^ is 
a stress on a stationary sensor, and K\ = ^irp s (d) 3 . 

(CP) 

Therefore, one expect that the distribution of er„ 



is approximately 2Ki 



Nv b 
-1-Af 



wider then the distribution 



of a 



(CV) 



This argument predicts that (width 



CP 



width cv )/width cv ~ 0.1 - 0.2, consistently with the re- 
sults shown in Fig. ll9f a') (in the estimate we use V = 10, 
based on the results shown in Fig.EJc)). We note that 
the above argument applies only to the normal stresses, 
since this is the direction of Vb- 



In principle, both normal and tangential stresses can 
also be affected by the changes in the volume fraction 
and in the velocity profiles due to the motion of the 
bottom wall. To estimate possible influence of this mo- 
tion, we note that it is characterized by the approxi- 
mate frequency of 1/(2At) = 457 Hz and the amplitude 
of VbAr/2 = 0.0005 (d). Despite high frequencies, the 
vibrations of such a small amplitude are not expected 
to significantly modify stress distributions, as also illus- 
trated also by the Fig. 1201 

Figure [2U| shows stress PDF's on typical size sensors 
using our typical averaging time At = T w /10. Com- 
paring to the stresses in CV simulations, Fig. I14f a1 and 
Fig. llBT a) . we do not see any significant differences. This 
is because At > At, therefore, during the averaging time 
the velocity of the bottom wall can change once or twice, 
resulting in "averaging out" the effect of the bottom wall 
motion. Therefore, we conclude that the model spec- 
ified by H28I30|I is effective in keeping constant stress, 
while not modifying significantly stress distributions on 
the time scales of interest. 



A. Bagnold scaling 

In this section we investigate the effect of the shearing 
velocity V on the stresses. The relation between average 
normal stress and V under certain conditions takes the 
form of Bagnold scaling 40 , where stresses increase with 
the square of the shear rate. This effect has been stud- 
ied in both CP and CV settings. 3 i 5 i 23 i 39 i 41 The quadratic 
dependence of the stress on the shearing rate is observed 
for the granular flows at high shearing rates or for lower 
volume fractions, i.e. in the systems where particles are 
not involved in multiple elastic deformations,- 

Figure I21f a) shows log-log plot of the mean normal 
stresses on the bottom wall versus V at CV setting for 
the system shown in Fig.^fa), i.e. with glued particles on 
the top wall and no oscillations. For configurations of v = 
40% and v — 50% these results confirm Bagnold scaling 
for the range of Vs between 1 rad/sec and 35 rad/sec. 
(Note that the v = 50% initial state was obtained from 
the typical configuration by adjusting the height of the 
cell according to ©). The data are fitted by the power 
function (a) = aV b with a = 0.015 and b = 2.20 ±0.2 for 
v = 40% and a = 0.201 and b = 2.18 ± 0.2 for v = 50%. 

Figure I21f b) shows scaled velocity profiles for different 
Vs in CV simulations. The profiles are higher for faster 
shearing. The origin of this dependence may be related 
to the effect of sidewalls on the sheared system. Indeed, 
when we shear the system with elastic and smooth side- 
walls, Fig. [f)| the profiles do not depend on V. Experi- 
mental studies^ also confirm the independence of veloc- 
ity profiles on V in the case boundary effects are signif- 
icantly minimized. To explain these results we consider 
the dependence of restitution parameters on the parti- 
cle's velocity in particle-sidewall interactions. The nor- 
mal coefficient of restitution Q depends on the normal 
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FIG. 19: (a) Stress on the base PDF for the case of CP simulations (circles), compared to CV results (triangles). Here 
At — At = T w /25 and A s = Abase for both sets of results, (b) The velocity of the bottom wall in the CP simulations recorded 
after some late time to in a steady state during 50 At intervals. 




FIG. 20: Constant pressure simulations: Stresses for the case without oscillations and with glued particles on the top wall. 
Size of sensors d s = 1.09(d). (a) (c) - sensors are on the top wall; (b) (d) - sensors are on the bottom wall; (a) (b) - normal 
stress PDF; (c) (d) - tangential stress PDF. 




FIG. 21: CV simulations, (a) Stress as a function of V in log-log plot. Circles: v = 40%. Triangles: v = 50%. Solid lines are 
the best fits of power function as explained in the text. The rest of results shown is calculated using v — 40%. (b) Velocity 
profiles; (c) Stress PDF, bottom wall, At = T w /25. Measurements are taken over last AT = 5000Ai = 5.47sec except for 
V = lrad/sec, where AT = 2sec; (d) Scaled energy plots for different Vs. 



velocity of colliding particle; however this velocity is not 
influenced significantly by V. On the other hand, the 
coefficient of tangential restitution, /3, in the case of slid- 
ing contacts |J5J depends stronger on V since it involves 
the ratio of normal to tangential velocity, and the later 
scales with V. As it can be seen from Q, lower ratio 
has the same effect on dissipative mechanics as lower co- 
efficient of friction, \i. Therefore, in the case of higher 
shear rates, the effect of sidewalls on the granular flow is 



reduced, explaining the velocity profiles in Fig. 121T b). 

The rest of the results shown in Fig. includes 
stresses and mean scaled energies. We discuss them later 
in this Section in the context of comparative study of 
constant volume and constant pressure boundary condi- 
tions. 

Next, we compare the effect of shearing velocity on 
the stresses for CV and CP configurations. In CP sim- 
ulations we choose the fixed value of the stress on the 
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FIG. 22: Constant pressure simulations. Position of the bottom wall, h, (in units of (d)) as a function of time for different 
shearing velocities. 



bottom wall given by l|29|l . Figure 12*2*1 shows the vertical 
position of the bottom wall h as a function of time for 
different V (h — corresponds to the distance from the 
top wall equal to 10.54(d)). As expected, we see that 
slower/faster shearing leads to smaller/larger stress, and 
to an adjustment of the wall position occurring over some 
time interval (which is influenced by the "spring" param- 
eters given by (|29I30|) 1. The stabilized (long time) wall 
positions are used to calculate v for each V Figure l^ST a) 
shows resulting z/s. We notice significant change in v as 

V is modified: v at high shearing (35 rad/sec) is about 
1.6 times smaller then v at low shearing (1 rad/sec). 

Figure I23f b) shows the effect of V on the velocity 
profiles. This result shows three types of different re- 
sponse to the increase of V in CP. For small Vs, V < 
2.5 rad/sec, the shearing is stronger for larger V. For in- 
termediate Vs, 2.5 rad/sec < V < 5 rad/sec the velocity 
profiles are almost the same, and for V > 5 rad/sec the 
shearing decreases with an increase of V These results 
can be understood in terms of the volume fraction, u, as 
a determining factor. For the smallest V (solid line in 
Fig. 123T b)) v is about 53%, therefore local jammed areas 
can be formed inside the sample reducing the overall mo- 
bility and shearing of the sample. These jammed areas 
can manifest themselves as "plateaus" in velocity profile. 
One such plateau can be seen in Fig. l23f b). solid line, be- 
tween y numbers 4 and 5. As v decreases with an increase 
of V, the sample gets more fluidized and the shearing im- 
proves. However, further increase of V and decrease of 
v makes the system more "compressible" : the dilation 
close to the top wall is stronger, reducing the momentum 
flux from the top wall to the bulk of the sample, resulting 
in weaker shearing. Therefore, the influence of V on a 
CP system is more complicated compared to a CV case, 
where velocity profiles are changing monotonously with 

V see Fig.ETTb). 

The distributions of the normal stress on the bottom 
wall are shown in Fig. l21f c) for the CV and in Fig. l23f c) 
for the CP. In both cases, CV and CP, the widths of 
distributions increase with a decrease of V. The distri- 
butions differ in the following: in CP case (Fig. l23t c'0 the 
distribution shows a double peak at low shearing around 



the prescribed value of the stress; This double peak is ab- 
sent in CV case. This double peak occurs since at low V 
the volume fraction is high, up to 52%. Thus the collision 
rate between the particles and the moving bottom wall 
is higher and, hence, larger momentum is transfered to 
the granular system. If this momentum is large enough, 
we have the situation in which the "restoring force" of 
the bottom wall brings the system in just one interval 
At from the overstressed state to under-stressed state or 
vice versa. Therefore, the velocity of wall, 128(1 , in steady 
state alternates its value between positive and negative 
without taking zero value, resulting in two peaks in stress 
distribution: one for overstressed state and another one 
for under-stressed state. 

Finally, Fig. I23f d') shows the time evolution of the 
mean scaled energy of the system. This figure shows that 
for large Vs, the equilibrating time decreases with an in- 
crease of V, while for low shearing this trend is opposite. 
This effect is due to the fact that the equilibrating time 
depends on the collision rate, which in turn is the func- 
tion of V and v. At high shearing v changes slowly with 
shearing, see Fig- I^ST a). so mostly V determines the col- 
lision rate, leading to the observed decrease of the equi- 
librating time with faster shearing. On the other hand, 
at slow shearing, small changes in V cause large changes 
in v. Here mostly v determines the collision rate, leading 
to the observed increase of the equilibrating time with 
faster shearing. 



VI. CONCLUSIONS 

The presented studies concentrate on 3D event driven 
simulations in the Couette geometry with top rotating 
wall and with physical boundary conditions in zero grav- 
ity. The velocity and volume fraction distributions are 
strongly dependent on various boundary conditions, such 
as: top shearing wall properties, side wall properties, 
presence of oscillations of the bottom wall, or intensity 
of shearing. The overall volume fraction of studied gran- 
ular system is 40 %, however, the shearing and vibrating 
walls impose the non-uniformity in volume fraction dis- 
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FIG. 23: CP simulations, (a) v as a function of V; (b) velocity profiles; (c) Stress PDF, bottom wall. These distributions are 
obtained from the stress data taken every At = At = T m /25 seconds over the last AT = 5000 At; (d) scaled energy plots. 



tribution with the formation of high-dense band (cluster) 
where local volume fraction reaches up to 60%. The clus- 
ter can respond to shearing in different ways. For exam- 
ple, while very rough and inelastic top wall, but without 
glued particles, cannot induce any significant shearing 
when side-walls are dissipating, it imposes a shear and 
high tangential velocities of all particles in the case of 
smooth and elastic side-walls. Typically, inside the clus- 
ter the shear rate is very small. However, the presence of 
the glued particles on the top wall considerably increases 
shear rate inside this cluster. 

The presence of oscillations of the bottom wall seems to 
have three effects: 1) the slippage velocity at the bottom 
increases. 2) the dense cluster is located further away 
from the bottom wall compared to the systems without 
oscillations. 3) equilibrating times to reach a steady state 
are shorter in the systems with oscillations due to the 
increased collision rate. 

The equilibrating dynamics involves attaining the 
steady state values of volume fractions and velocities. 
The time scales of these two equilibrating processes are 
different: in most cases, steady state profile in volume 
fraction is attained very fast, while velocity profiles reach 
their final shape much later. However, the volume frac- 
tion distribution takes longer to reach a steady state 
when the rate of energy input due to shearing is low. This 
occurs in the systems without glued particles and without 
oscillations. For example, in so called "delayed dynam- 
ics" regime, described in Sec. IHI CI the time needed for 
the cluster to accumulate enough energy to dis-attach 
from the bottom wall can be very long. However, once 
this happens, the velocity profile changes drastically. 
More generally, these simulations show that the velocities 
and volume fractions strongly depend on boundary con- 
ditions: simulations using e.g. periodic boundary condi- 
tions may miss a number of interesting effects presented 
in this work. 



We have also analyzed stress and stress fluctuations on 
the boundaries. We find that the distribution of normal 
stresses in a system of 40% volume fraction and rapid 
granular flow can be described by the same functional 
form (|I2|I as used for stress distribution in static dense 
system, i.e. with power-law increase for small stresses 
and exponential decrease for large ones. Considering this 
functional form as a useful guide, we predict the behavior 
of main characteristics of stress distributions as a func- 
tion of sensor area and averaging time. These predic- 
tions describe well the width of distributions for small 
sensors and short At's. However, they fail for large sen- 
sors and long averaging times, suggesting the existence of 
additional correlations not accounted by our functional 
form. The correlations may be related to the existence 
of "denser structures" in the simulated systems. 

We simulate "constant pressure" boundary condition 
and we discuss the differences between the results in con- 
stant volume and constant pressure settings. A key ob- 
servation here is very different reaction to an increase of 
shearing velocity in the system with constant pressure 
boundary condition compared to the systems with con- 
stant volume boundary condition. In the first case the 
shearing velocities are increasing because of the veloc- 
ity dependence of the coefficients of restitution while in 
the second case they are initially increasing and then de- 
creasing because of the changes of volume fraction with 
imposed shear. 
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